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Abstract 

The bifurcation structure of coupled periodically driven double-well Duffing 
oscillators is investigated as a function of the strength of the driving force 
/ and its frequency O. We first examine the stability of the steady state in 
linear response, and classify the different types of bifurcation likely to occur 
in this model. We then explore the complex behaviour associated with these 
bifurcations numerically. Our results show many striking departures from the 
behaviour of coupled driven Duffing Oscillators with single well-potentials, as 
characterised by Kozlowski et al [1]. In addition to the well known routes to 
chaos already encountered in a one-dimensional Duffing oscillator, our model 
exhibits imbricated period-doubling of both types, symmetry-breaking, sud- 
den chaos and a great abundance of Hopf bifurcations, many of which occur 
more than once for a given driving frequency. We explore the chaotic be- 
haviour of our model using two indicators, namely Lyapunov exponents and 
the power spectrum. Poincare cross-sections and phase portraits are also 
plotted to show the manifestation of coexisting periodic and chaotic attrac- 
tors including the destruction of T 2 tori doubling. 
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I. INTRODUCTION 



In recent years, a large number of theoretical calculations, numerical simulations and 
experiments has been carried out on systems of coupled anharmonic oscillators which provide 
fundamental models of the dynamical problems in several disciplines. Such models describe 
a wide range of processes, helping in general to understand the routes to chaos that take 
place in biological, chemical, physical, electromechanical and electronical systems. Here we 
refer only to recent papers close to the subject of our own[l-18]. Coupling may arise between 
oscillators of the same types or oscillators of different types. Amongst these systems, the 
most intensively investigated examples are the Duffing oscillators [1-8] (strictly dissipative 
systems which converge to quiescence when not driven), the Van der Pol oscillators [9- 11] 
(self-excited systems which final behavior when not driven is a limit cycle) and the coupled 
Van der Pol - Duffing oscillators [12-13]. 

In this paper we present the results of an investigation of two identical coupled double-well 
Duffing oscillators subjected to a periodically driven force. Our study was stimulated by 
the earlier work of Kozlowski et al. [1] who have studied the model of the same type, but 
considered only single-well potentials. Using bifurcation diagrams and phase diagrams these 
authors have demonstrated that the global pattern of bifurcation curves in parameter space 
consists of repeated subpatterns similar to the superstructure found in a one-dimensional 
driven Duffing oscillator. They have also proven the existence of a Hopf Bifurcation which 
does not appear in a model of one-dimensional Duffing oscillator. Other interesting results 
have been reported in the context of coupled Duffing Oscillators. In ref. [2] Kunick et 
al showed that two coupled oscillators can behave regularly with nonzero coupling and 
chaotically with zero coupling. In a model of two coupled Duffing oscillators driven with 
incommensurate frequencies and coupled additively, Stagliano et al [3] observed the period- 
doubling of an attracting T 2 torus and its destruction in parameter space. Migration control 
in two Duffing oscillators by open-plus-closed-loop control method and adaptive control 
algorithm were studied by Paul et al [4]. These authors, in another system of two coupled 



2 



Duffing oscillators, have also investigated a basin of attraction with several coexisting chaotic 
attractors and synchronization of chaos [5]. Yagasaki [6], after having described an Auto 
driver Hommap for numerical analysis of homoclinic in maps of periodically forced Duffing 
systems, modified a version of the homoclinic Melnikov method for orbits homoclinic to two 
types of periodic orbits. These theories have been successfully applied to a weakly coupled 
Duffing oscillators. Moreover Yin et al [7] investigated the effect of phase difference in the 
driving forces of two coupled identical Duffing oscillators. They observed desynchronization 
and lag synchronization of chaotic attractors. Bifurcation behaviors, showing a kind of Hopf 
bifurcation, from synchronous chaos of a chain of Duffing oscillators have been explored by 
Ma Wen-Qi et al. [8] using the generalized winding number in tangent space. 
The equations of motion for the two identical coupled double-well Duffing oscillators that 
we are interested in are the following: 



where a is the damping parameter, k the coupling parameter, / and Vt the amplitude and 
the frequency of the driving force, respectively. In system (1), we represent the oscillators 
by the state variables (x,dx/dt) as oscillator A or subsystem A subjected to the periodic 
force f(t) = /cos(fii), which is coupled to the other one with the state variables (y,dy/dt) 
as oscillator B or subsystem B. Here / and are the control parameters of our system. 
As evident, for k = 0, equation (1) describes two uncoupled Duffing Oscillators with their 
motion governed by a and /. The coupling considered can be interpreted as a perturbation 
of each oscillator through a signal proportional to the difference of their amplitudes. A 
natural question to ask is how chaotic attractors arise as these parameters are varied. In the 
case of one dimensional double-well Duffing oscillator, quasiperiodic and period-doubling 
routes to chaos have been found [19-27]. 

The rest of the paper is structured as follows. In section II, we establish a linear stability 
analysis of the system projected on the Poincare map. Section III is devoted to the survey 



d 2 x dx 



x 3 + k(y - x) + f cos(fit) 



(1) 




y 3 - k(y - x ) 
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of bifurcation diagrams and the chaotic behavior is characterized in section IV. We end with 
the conclusions of the work in section V. 



II. STABILITY ANALYSIS 

For stability analysis and even for purposes of numerical simulations, it is convenient 
to transform the second order differential equation (1) into an autonomous system of first 
order differential equations of the following form: 

dx 

— — v 
dt 

dv 

di 
dy 

dt 
dw 

~dt 

de a 



= —av + x — x 3 + k(y — x) + / cos(0) 

(2) 

= —aw + y — y 3 — k(y — x) 



dt 2tt 



or equivalently 



where T = 2n/Vt is the period of the driving force, 6 the cyclic variable, V(x, v, y, w, 9) an 
autonomous vector field and //(a, k, /, Q) an element of the parameter space. This system 
generates a flow = {4> T } on the phase space R 4 x S 1 and there exists a global map: 

P : S c > S c 

V p {x,v,y,w) -> P(V P ) = {0 T }| Sc (x, w a.) 

with 6q being a constant determining the location of the Poincare cross-section and 
(x, v, y, w) the coordinates of the attractors in the Poincare cross-section J2 C defined by: 
X) c — v i Di w i 8) G R 4 x S 1 ]^ = ^ }- System (2) is symmetric since the transformation 
S : (x,v,y,w,9) — > (— x, —v, —y, —w,9 + n) leaves (2) invariant. This system can support 
symmetric orbits and also asymmetric ones which are not invariant with the transformation 



S. We use a perturbation analysis to proceed analysis of solution stability in the poincare 
map £ c . This method consists of differentiation and translation. Thus if we add a small 
perturbation 8V p (5x, 5v, 5y, 5w) to the steady state V p o(xo,vo,yo,wo) we obtain equations 
for perturbed motion 



x = x + 5x, y = yo + Sy, v = v + 8v, w = w + 5w 



(3) 



Substituting equations (3) into equations (2) and developing into powers of the above small 
perturbations yields the following variational differential equation: 



dSV p 
dt 



DF(V p0 )SV p , 



(4) 



with 



DF(V p0 ) 



1 — k — 3x1 —a 








k 







1 



V 



/ 





k l-k-3y% -a 

where DF(V p0 ) is the Jacobian 4x4 matrix describing the vector field along the solution 
SV p (t) and VpQ being an equilibrium point or the steady state. Next for the very small 
amplitudes, we neglect in the matrix DF(V p o) the powers higher than one, that is, Xq and 
?/q. The solution after one period T of the oscillations in the linearised Poincare map is 
simply expressed as: 



SV P (T) = SVpo exp(DF(V p0 )T), SV p (0) = 5V p0 



(5) 



where DF(V p o) is the monodromy matrix of a periodic orbit connecting arbitrary infinites- 
imal variations in the initial conditions 8V p0 with corresponding changes 5V P (T), after one 
period T. The stability of the periodic motion is therefore determined according to the 
real parts of the roots of the characteristic equation det(DF(V p o) — IX) = written in the 
following form: 



A 4 + A 3 \ 3 + A 2 \ 2 + A 1 \ + A = 



(6) 



where (Ai),i = 0,1,2,3 are the coefficients depending on the two parameters k, a and A 
= (Aj) the eigenvalues of DF(V p0 ). The roots of equation (6) are obtained using the Bairstow- 
Newton-Raphson algorithm [28] and for different values of k G [—5,5] and a G [—0.2,0.2]. 
The results as the spectrum of the eigenvalues are depicted in a complex plane C (see fig.l). 
This figure enables us to get an idea of different types of bifurcation likely to appear in the 
system. Since DF(V p o) is a real matrix, complex eigenvalues occur in a complex conjugate 
pairs responsible of the symmetry observed along the real axis. Thus if A, is real, it is 
clear from equation (5) that the eigenvalues is the rate of contraction (Aj < 0) or expansion 
(Aj > 0) near the steady state. Next if Aj is complex, the real part of Aj gives the rate 
of contraction (ite(Aj) < 0) or expansion (ite(Aj) > 0) of the spiral while the imaginary 
part Im(\i) is of the frequency rotation. This can be well understood with the following 
expression of the eigenvalues of the linearized Poincare map: 

<7j = exp(Re(\i))[cos(Im(Xi)T) + jsin{Im{\)T)] (7) 

Globally, if Re{\) < for all Aj, then all sufficiently small perturbations tend toward zero 
as t goes to infinity and the steady state (nodes (n), saddle nodes (sn), spiral (sp)) is stable. 
If Re{\) > for all Aj, then any small perturbation grows with the time and the steady 
state (n, sn, sp) is unstable (A stable or unstable equilibrium state with no complex eigen- 
values is often called node). On the other hand if there exists % and I such that Re(Xi) < 
and Re(\i) > 0, therefore the equilibrium state is non-stable (A non-stable equilibrium 
state is often called saddle; an equilibrium point whose eigenvalues all have a non zero real 
part are called hyperbolic [29]). Having regards to the precedent with again a look on the 
spectrum of fig.l, it follows that our system of the given parameter k and a can undergo 
many types of bifurcations namely: saddle-node (ite(Aj) = 1), period-doubling (Aj = — 1), 
Hopf bifurcations (Aj = 7 ± j(5 , (j 2 = —1) with 7 < and the motion has the form 
of exponential increasing beats with period 2n/j3 [17]) and of course symmetry-breaking 
bifurcation which is often a prerequisite for the first period-doubling bifurcation [30]. Ex- 
cept Hopf bifurcation, such bifurcations have been successfully found in one-dimensional 
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double- well Duffing oscillators subjected to a periodically driven force [19,21,26] using the 
eigenvalues called Floquet multipliers of the linearized Poincare map. In a coupled identical 
single-well Duffing oscillators, Kozlowski et al [1] showed by linear analysis a similar eigen- 
values scenario and they conjectured that the resulting alternating bifurcation sequences 
have to be expected for all systems of coupled strictly dissipative oscillators that are driven 
periodically. The above analysis allows us to know how the steady state becomes locally 
unstable and to be aware of the type of bifurcation expected in the system. Here, no explod- 
ing amplitude is possible since the shape of the potential offers globally bounded solutions 
(V(x,y) = -(l/2)x 2 - (l/2)y 2 ) + (l/4)x 4 + (l/% 4 ) -> +00 as \x\,\y\^oo). 

III. LOCAL BIFURCATIONS 

The aim of this section is to seek numerically the routes which lead to chaotic solutions 
of the system when the control parameter Q evolves, for different values of the driving force 
/. Bifurcation diagrams are very good for numerical as well as for experimental studies 
when there is a tunable parameter. When a control parameter is varied and bifurcation 
takes place a qualitative change of the system happens. For the numerical computations 
of these diagrams the control parameter Q is increased from an initial value f2j to a final 
value £lf and then decreased from to ^ in a very small step (Af2 = 1CT 5 ). The last 
computed cyclic point for a given value of Q is always used as a new initial value for the 
next value of Q. Starting with initial conditions (x,v) = (1,0), (y, w) = (1,0) at fij, system 
(2) is integrated, using the standard Runge Kutta Algorithm [28], for 100 periods of the 
driving force until the transient has died out; the trajectories expected are attractors and 
the local calculation error is sufficiently small. Then to find out whether the trajectory 
is periodic (quasiperiodic) or chaotic, the system is integrated for the next 180 periods in 
order to catch the maximum of coexisting attractors. This procedure allows us to make sure 
that coexisting attractors are realised due to different initial conditions which are handed 
throughout from one parameter to another. We chose k = 5 and a — 0.1 for all examples 
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in the subsequent section. The bifurcation diagrams obtained here show the projection of 
the attractors in the Poincare section onto x or v versus the control parameter Q and for 
different values of /. In a one-dimensional double- well Duffing oscillator [19-21,23-27], it 
is shown that chaotic solutions arise through quasiperiodic, period-doubling cascade with a 
sequence of symmetry-breaking and saddle-nodes. Hopf bifurcations have not been found, 
but are abundant in our model. 

We computed a great number of diagrams, and some of them showing the main features 
of the system are presented in figures 2-4. It is important to note that no hysteresis has 
been detected while checking all these diagrams. We can therefore avoid any confusion 
over the periods of attractors when reading this type of diagram; for instance one single 
period-2 yields to 2 points for one frequency value whereas two coexisting attractors each of 
period- 1 also yield to two points for one frequency. These diagrams show a great number of 
coexisting attractors (chaotic domains) intermingled with imbricated windows made up of 
periodic solutions of different periodicity, period-doubling of both types, sudden Chaos and 
Hopf bifurcation. The so called chaotic domains (or chaotic sea) will be characterized in the 
following section. 

At lowest frequencies (0 < f2 < 0.4) quasiperiodic as well as chaotic solutions occur together 
with a number of symmetry-breaking (sb), period doubling (pd) and saddle-nodes (sn). In 
fig.2(a) where / = 15, it appears in the neighborhood of f2 = 0.385, 6 pairs of sb followed 
by 3 reversed pd till the chaotic domain. When the driving force is increased / = 20, 
overlapping pd cascade and sn appear as can be seen in fig.2(b). Similar features as in the 
latter case are also observed for small values of / (/ < 15). In other ranges of the driving 
frequency bigger than Q = 0.4 the aformentioned features still exist and additional ones 
such as resonances and Hopf bifurcations appear in the system; the results are displayed in 
figures 3 and 4. In figure 3 it can be seen throughout overlapping pd cascade of both types, 
sb, resonances (R) and even sudden chaos (SC). When / = 4.5 in fig. 3(f) at Q ~ 0.785, 
period-8 attractor undergoes reversed period-doubling cascade and period-2 is created as 
from f2 ~ 0.805. Besides, for / > 7 appear resonances (R) in fig.3(a,b,d,e), fig.6(a) and 
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sudden chaos in fig. 3(c). However in fig. 4 when / becomes much larger and Q > 0.5, 
Hopf Bifurcation (H) is abundant and often appear more than once at a given frequency. 
This notation —nH— means that there are n, (n < 3), integer Hopf bifurcation at a given 
frequency. One can therefore observe -H- in (a), -3H- in (b), -2H-..-2H- in (c,d), -2H-..-3H-..- 
2H- in (e,f) and -3H-3H- in (g,h). In fig. 4(h) for instance, at f2 = 0.5925, a period-6 attractor 
undergoes reversed period-doubling to become a period-3; next after 3 pairs of symmetry- 
breaking near £7 ~ 0.6, appears again period-3 attractor as from ~ 0.601. And then at 
f2 = 0.603 this attractor of period-3 undergoes 3 Hopf bifurcations. Windows of periodic 
solutions separated by quasiperiodic ones are clearly visible. The Poincare cross-section of 
quasiperiodic attractors of this kind (-3H-) consists of three invariant tori (destroyed or not). 
Thus at any f2 where —nH— exist, the projection of the attractor corresponds to n Tori 
on the Poincare cross-section. The torus attractor arises from Hopf bifurcation which is 
a bifurcation from a fixed point to an invariant curve. At the transition to quasiperiodic 
motion, this fixed point (projection of the limit cycle in to the cross-section) loses its stability 
and gives birth to an invariant circle, which is a cross-section of the torus in the flow. The 
three pairs of sb appearing at the neigborhood of ~ 0.6 (fig.4(g,h)) seem to be apparently 
a window of period-6 attractor. In fact it is just a period-3 attractor as can be illustrated 
further in fig.8(b).This orbit, born asymmetrically, coexists with its inversion symmetric 
orbit of the same period-3 and are caught simultaneously. 

Next we pursue our investigation for large values of f2 (Q > 3) and for different values of 
/. It turns out that the features observed so far have disappeared and the system becomes 
regular providing attractors with a small number of periods as can be seen in figure 5. The 
phase diagrams which often show bifurcations, curves, surfaces and which is out of the scope 
of this paper could also help to match the features obtained in this analysis. 
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IV. CHAOTIC BEHAVIOR 



To better support the results obtained above, we intend in this section to characterize 
chaotic behavior. For this purpose some indicators are used namely the Lyapunov exponents 
and the power spectrum density. Also, Poincare cross-sections and phase portraits are shown 
to illustrate some attractors present in the system. The Lyapunov exponent is one of the 
most important tools for understanding chaotic behavior, obtained by examining "a very 
sensitive dependence of initial conditions" . In particular it is generally well established that 
a rigorous measure of chaos may be given in terms of Lyapunov spectrum of the dynamical 
system. A positive Lyapunov exponent is characteric of chaos while zero and negative 
values of the exponent signify a marginally stable or quasiperiodic orbit and periodic orbit, 
respectively. Solving numerically the variational equation (4) together with system (2) 
we calculate the maximum Lyapunov exponent \ m ax i n the Poincare cross-section with 
^max = lim T ^+oo 7 hi(||L(r) ||), where ||£(r)|| = (5x 2 + 5v 2 + by 2 + 5w 2 ^j ^ . The spectrum 
is displayed in fig. 6 together with a bifurcation diagram. Windows of chaotic domains 
and periodic solutions are clearly justified. Power Spectrum density using Fast Fourier 
Transformation (FFT) has also been computed with some values from fig. 6 for a periodic 
solution (Q = 0.495, / = 7) and a chaotic solution (Q = 0.47,/ = 7) (see fig.7). Fig.8 
shows the phase portraits and Poincare cross-section. In fig.8 (a,b) plotted are the phase 
portraits, related to fig.4(g,h), showing a period-3 attractor for / = 30, Q = 0.6025(a) 
and a coexistence of two asymmetric attractors, each of period-3 for / = 30, = 0.6 (b). 
The dots (•) represent the corresponding points in the Poincare cross-sections. In fig.8 (c,d) 
two chaotic attractors are plotted in the Poincare cross-section for Q = 0.5, / = 2 and 
for f2 = 0.375, / = 15, respectively. Furthermore, we observe in fig.8 (e,f) three T 2 tori 
doubling at f2 = 0.605, / = 30 and at Q = 1.14, / = 15, two T 2 tori doubling (destroyed) 
corresponding for -3H- and -2H-, respectively (values from fig.4(c,h)). In the torus route to 
chaos, the original torus appears to split into two circles at the torus doubling bifurcation 
point. This route is reminiscent to the period-doubling route to chaos, although there are 
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finite number of torus doublings before the onset of chaos motion. Also, the torus doubling 
route to chaos is a higher-dimensional phenomenon, requiring at least a four- dimensional 
flow, or a three-dimensional map. It is not observed in one-dimensional maps, unlike the 
period-doubling route to chaos [30] . 

V. CONCLUSION 

To sum up we have investigated coupled double-well Duffing oscillators subjected to a 
periodic driving force. Our main conclusion is that the shape of the potential in the coupled 
Duffing Oscillators system has a profound influence on the routes to chaos taken by the sys- 
tem. By a linear analysis we have examined the stability of the steady state of the system 
leading to different types of bifurcation likely to appear for k e [—5,5] and a G [—0.2,0.2]. 
With the help of bifurcation diagrams, it has been shown that the bifurcation structure, 
which is complicated, depends strongly on the values of the control parameter Q. For a 
simple electronic realization of our model, it is obvious that the oscillators may be very 
difficult to control for small values of f2 due to the large windows for chaos. This system 
becomes regular for large values of Q (Q > 3). Apart from the routes already encountered in 
a one-dimensional double-well Duffing oscillator such as quasiperiodic and period-doubling 
cascade, sudden chaos and mostly Hopf bifurcations have been newly found in this work. 
Besides, resonances (R), imbricated sb, sn, pd intermingled by chaotic domains are ob- 
served, rending the structure highly chaotic. Two and three T 2 tori doubling, have been 
also observed. To characterize chaotic behavior of this system, Lyapunov exponents and the 
power spectrum density using FFT have been employed. Poincare cross-sections and phase 
portraits have been also used for better visualisation of periodic and chaotic attractors. An 
analytical approach using approximative technics could be an interesting topic for future 
work, and would probably shed further light on the rich behaviour of this model. 
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VII. FIGURE CAPTION 



Figure 1: Eigenvalues spectrum for the stability of the system with a G [—0.2, 0.2] and 
k E [—5,5]. Note that when a < 0, the eigenvalues are in the half-plane left and in 
the half-plane right otherwise. 

Figure 2: Bifurcation diagrams showing the first coordinate x for the driving amplitude 
/ = 15 (a) and the second coordinate v for / = 20 (b) of the Poincare section versus the 
driving frequency ft. Windows of periodic solutions and chaotic domains are visible. 
Near ft = 0.385 (a), 6 pairs of sb occur followed by the reversed pd till the chaotic 
domain. Imbricated pd, sb, sn are mentioned. Note that similar behavior as (b) is 
observed when ft < 0.4 and for small values of / (/ < 15). 

Figure 3: Bifurcation diagrams for the driving force / = 4.5 (f), / = 15 (a,b,e), 
/ = 20 (d), / = 25 (c) with ft > 0.4. Throughout imbricated sb, sn, pd of both types 
and chaotic domains are observed. Besides, resonances (R) appear at ft ~ 0.475 (a,b), 
ft ~ 0.445 (d), ft ~ 0.555 (e)and sudden chaos at ft ~ 0.408 (c). 

Figure 4: Bifurcation diagrams for / = 30 (a,g,h), / = 25 (b), / = 15 (c,d), / = 20 (e,f) 
showing many Hopf bifurcations (H). One can observe -H-(a), -3H-(b), -2H-..-2H (c,d),- 
2H-..-3H-..-2H-(e,f), -3H-3H- (g,h) (see also text) 

Figure 5: Bifurcation diagram for / = 5 showing attractors of small number of periods. 
Note that the system becomes regular and this behavior is observed for ft > 3 and for 
any value of /. 

Figure 6: Bifurcation diagram for / = 15 (a) and the corresponding Lyapunov spec- 
trum (b) in the Poincare map. windows of chaotic domains are clearly justified with 
positive values of \ max - 

Figure 7: Power spectra densities (PSD). The values of / = 7 and ft are from the 
bifurcation diagrams shown in fig.6 with ft = 0.495 (a) periodic motion and ft = 

15 



0.47 (b) chaotic motion. 

• Figure 8 : Phase portraits and projection of the attractors in the Poincare section onto 
the two first coordinates x,v of the Poincare cross- sect ion. At Q, = 0.6025 (a), one 
period-3 orbit and at Jl = 0.6 (b) coexistence of two asymmetric orbits of the same 
period-3. The values of Q are from fig. 4 (g,h) and / = 30. The dots (•) represent the 
corresponding points in the Poincare cross-section. Two chaotic attractors shown in 
(c,d) for Q = 0.5, / = 2 and Q = 0.375, / = 15, respectively. Three T 2 tori doubling 
for n = 0.605, / = 30 (e) and two T 2 tori doubling (destroyed) for SI = 1.14, / = 15 
(f). The values of Q are from fig.4(c,h). 
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